Methods Overview

  • Multidimensional IRT used to co-calibrate HRS and HCAP cognitive measures
    • Merged HRS-HCAP data set from 2016 (HRS wave 13) used for co-calibration
  • Measures evaluated
    • TICS3 (wave 3)
      • Word List Recall (immediate and delayed),
      • Count Backwards
      • Serial 7’s
      • Orientation to Time
      • Naming
    • TICS10 (wave 3 + Number Series)
    • TICS3 + Number Series
    • TICS13 (in TICS variants) / TICS (Comparisons with non-TICS measures)
    • TICS3 + Number Series + Animal Fluency
    • TICS16
    • TICS13 + Trail Making A & B
    • MMSE
    • HCAP
    • MMSE+TICS+HCAP
  • Simulation of cognitive trajectories for different measures
    • True cognition simulated based on empirical modeling of cognitive trajectories from a large sample of demographically and racially diverse older adults
      • UC Davis ADRC
      • Kaiser cohort studies
        • KHANDLE
        • STAR
        • Life After 90
    • Simulated item response datasets created using true cognition values and item parameters from multidimensional IRT co-calibration

Data

# ********************************** Data **************************************


hcog <- read.csv("~/Psychometrics Conference/2024/Simulation WG/Data/HRS_Merge3.csv")
names(hcog) <- tolower(names(hcog))


# MMSE

hcog <- hcog %>% mutate(
    ornttime_mmse = r1orient_time,
    orntplace_mmse = r1orient_place,
    imrec_mmse = r1imrc3,
    delrec_mmse = r1dlrc3,
    spellbck_mmse = backspel,
    naming_mmse = r1object,
    phrase_mmse = case_when(
        !is.na(r1repeat) ~ case_when(
            r1repeat == 1 ~ 1,
            r1repeat == 5 ~ 0
        )
    ),
    commfoll_mmse = case_when(
        !is.na(r1combfol) ~ case_when(
            r1combfol %in% c(1,2) ~ 1,
            r1combfol == 5 ~ 0
        )
    ),
    writesent_mmse = case_when(
        !is.na(r1senten) ~ case_when(
            r1senten %in% c(1,2) ~ 1,
            r1senten == 5 ~ 0
        )
    ),
    draw_mmse = case_when(
        !is.na(r1draw) ~ case_when(
            r1draw == 1 ~ 1,
            r1draw == 5 ~ 0
        )
    ),
    comm3step_mmse = follow3_step,
    mmsetot_mmse = r1mmse_score
)

# HRS TICS

hcog <- hcog %>% mutate(
    wlimrc_tics = pd174,
    wldlrc_tics = pd184,
    # count_tics = pd170,
    animfl_tics = pd196,
    numser_tics = pd216,
    serial7_tics = ser7sum,
    ornttime_tics = case_when(
        !is.na(pd151) ~ case_when(
            pd151 == 1 ~ 1,
            pd151 == 5 ~ 0
        )
    ) +
        case_when(
            !is.na(pd152) ~ case_when(
                pd152 == 1 ~ 1,
                pd152 == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(pd153) ~ case_when(
                pd153 == 1 ~ 1,
                pd153 == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(pd154) ~ case_when(
                pd154 == 1 ~ 1,
                pd154 == 5 ~ 0
            )
        ),
    naming_tics = 
        case_when(
            !is.na(pd155) ~ case_when(
                pd155 == 1 ~ 1,
                pd155 == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(pd156) ~ case_when(
                pd156 == 1 ~ 1,
                pd156 == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(pd157) ~ case_when(
                pd157 == 1 ~ 1,
                pd157 == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(pd158) ~ case_when(
                pd158 == 1 ~ 1,
                pd158 == 5 ~ 0
            )
        ),
    cntbck_tics = pd170 - ornttime_tics - naming_tics
)

# HCAP

hcog <- hcog %>% mutate(
    naming_hcap = 
        case_when(
            !is.na(r1scis) ~ case_when(
                r1scis == 1 ~ 1,
                r1scis == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(r1cactus) ~ case_when(
                r1cactus == 1 ~ 1,
                r1cactus == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(r1pres) ~ case_when(
                r1pres == 1 ~ 1,
                r1pres == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(r1elbow) ~ case_when(
                r1elbow == 1 ~ 1,
                r1elbow == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(r1hammer) ~ case_when(
                r1hammer == 1 ~ 1,
                r1hammer == 5 ~ 0
            )
        ),
    ceradtot_hcap = r1word_total,
    ceraddel_hcap = r1word_dscore,
    ceradrcg_hcap = r1wlrec_totscore,
    animfl_hcap = r1verbal_score,
    bmim_hcap = r1bm_immscore,
    bmdel_hcap = r1bm_delscore,
    lmim_hcap = r1lmb_immscore,
    lmdel_hcap = r1lmb_delscore,
    lmrcg_hcap = r1lmb_recoscore,
    conprxim_hcap = r1cp_score,
    conprxdel_hcap = r1cpdel_score,
    dsmt_hcap = r1dig_score,
    numser_hcap = r1ns_score,
    ravens_hcap = r1rv_score,
    trailsa_hcap = r1tma_score,
    trailsb_hcap = r1tmb_score,
    spatial_hcap =
        case_when(
            !is.na(r1store) ~ case_when(
                r1store == 1 ~ 1,
                r1store == 5 ~ 0
            )
        ) +
        case_when(
            !is.na(r1point) ~ case_when(
                r1point == 1 ~ 1,
                r1point == 5 ~ 0
            )
        )
)

mmse_nms <- c("ornttime_mmse","orntplace_mmse","imrec_mmse","delrec_mmse",
    "spellbck_mmse","naming_mmse","phrase_mmse","commfoll_mmse","writesent_mmse",
    "draw_mmse","comm3step_mmse")

tics_nms <- c("wlimrc_tics","wldlrc_tics","animfl_tics",
    "numser_tics","serial7_tics","ornttime_tics","naming_tics","cntbck_tics")

hcap_nms <- c("naming_hcap","ceradtot_hcap","ceraddel_hcap","ceradrcg_hcap",
    "animfl_hcap","bmim_hcap","bmdel_hcap","lmim_hcap","lmdel_hcap","lmrcg_hcap",
    "conprxim_hcap","conprxdel_hcap","dsmt_hcap","numser_hcap","ravens_hcap",
    "trailsa_hcap","trailsb_hcap","spatial_hcap")

### Recode MMSE

# summary(hcog[,mmse_nms])

hcog$imrec_mmse <- missCode(hcog,"imrec_mmse",0,3)
hcog$delrec_mmse <- missCode(hcog,"delrec_mmse",0,3)

# table(hcog$ornttime_mmse)
# table(hcog$orntplace_mmse)
# table(hcog$imrec_mmse)
# table(hcog$delrec_mmse)
# table(hcog$spellbck_mmse)
# table(hcog$naming_mmse)
# table(hcog$phrase_mmse)
# table(hcog$commfoll_mmse)
# table(hcog$writesent_mmse)
# table(hcog$draw_mmse)
# table(hcog$comm3step_mmse)

hcog <- hcog %>% mutate(
    ornttime_mmser = ornttime_mmse,
    orntplace_mmser = orntplace_mmse,
    imrec_mmser = case_when(
        imrec_mmse %in% c(0,1) ~ 0,
        TRUE ~ imrec_mmse - 1
    ),
    delrec_mmser = delrec_mmse,
    spellbck_mmser = spellbck_mmse,
    naming_mmser = naming_mmse,
    phrase_mmser = phrase_mmse,
    commfoll_mmser = commfoll_mmse,
    writesent_mmser = writesent_mmse,
    draw_mmser = draw_mmse,
    comm3step_mmser = comm3step_mmse
)

# table(hcog$imrec_mmse,hcog$imrec_mmser)

### End recode MMSE

### Recode TICS

# summary(hcog[,tics_nms])
# 
# table(hcog$wlimrc_tics)
# table(hcog$wldlrc_tics)
# table(hcog$cntbck_tics)
# table(hcog$count_tics)
# table(hcog$animfl_tics)
# table(hcog$numser_tics)
# table(hcog$serial7_tics)
# table(hcog$ornttime_tics)
# table(hcog$naming_tics)

hcog <- hcog %>% mutate(
    wlimrc_ticsr = case_when(
        wlimrc_tics %in% c(9,10) ~ 9,
        TRUE ~ wlimrc_tics
    ),
    wldlrc_ticsr = case_when(
        wldlrc_tics %in% c(9,10) ~ 9,
        TRUE ~ wldlrc_tics
    ),
    cntbck_ticsr = case_when(
        cntbck_tics %in% c(0,1) ~ 0,
        TRUE ~ cntbck_tics - 1
    ),
    # count_ticsr = case_when(
    #     count_tics %in% c(0,1,2) ~ 0,
    #     TRUE ~ count_tics - 2
    # ),
    ornttime_ticsr = case_when(
        ornttime_tics %in% c(0,1) ~ 0,
        TRUE ~ ornttime_tics - 1
    ),
    naming_ticsr = case_when(
        naming_tics %in% c(0,1,2) ~ 0,
        TRUE ~ naming_tics - 2
    ),
    numser_ticsr = numser_tics,
    serial7_ticsr = serial7_tics
)

# table(hcog$wlimrc_tics,hcog$wlimrc_ticsr)
# table(hcog$wldlrc_tics,hcog$wldlrc_ticsr)
# table(hcog$count_tics,hcog$count_ticsr)
# table(hcog$ornttime_tics,hcog$ornttime_ticsr)
# table(hcog$naming_tics,hcog$naming_ticsr)

hcog <- recodeOrdinal(df = hcog, varlist_orig = "animfl_tics", varlist_tr = "animfl_ticsr")

# table(hcog$animfl_tics,hcog$animfl_ticsr)

### End recode TICS


### Recode HCAP

# summary(hcog[,hcap_nms])

hcog$trailsa_hcap <- missCode(hcog,"trailsa_hcap",0,300)
hcog$trailsb_hcap <- missCode(hcog,"trailsb_hcap",0,300)

hcog$trailsa_hcapi <- -hcog$trailsa_hcap
hcog$trailsb_hcapi <- -hcog$trailsb_hcap

# summary(hcog[,c("trailsa_hcapi","trailsb_hcapi")])
# 
# table(hcog$naming_hcap)
# table(hcog$ceradtot_hcap)
# table(hcog$ceraddel_hcap)
# table(hcog$ceradrcg_hcap)
# table(hcog$animfl_hcap)
# table(hcog$bmim_hcap)
# table(hcog$bmdel_hcap)
# table(hcog$lmim_hcap)
# table(hcog$lmdel_hcap)
# table(hcog$lmrcg_hcap)
# table(hcog$conprxim_hcap)
# table(hcog$conprxdel_hcap)
# table(hcog$dsmt_hcap)
# table(hcog$numser_hcap)
# table(hcog$ravens_hcap)
# table(hcog$trailsa_hcap)
# table(hcog$trailsb_hcap)
# table(hcog$spatial_hcap)

hcog <- hcog %>% mutate(
    naming_hcapr = case_when(
        naming_hcap %in% c(0,1,2,3) ~ 1,
        TRUE ~ naming_hcap - 2
    ),
    spatial_hcapr = spatial_hcap + 1
)

varlist_orig <- c("ceradtot_hcap","ceraddel_hcap","ceradrcg_hcap",
    "animfl_hcap","bmim_hcap","bmdel_hcap","lmim_hcap","lmdel_hcap","lmrcg_hcap",
    "conprxim_hcap","conprxdel_hcap","dsmt_hcap","numser_hcap","ravens_hcap",
    "trailsa_hcapi","trailsb_hcapi")

varlist_tr <- paste0(varlist_orig,"r")

hcog <- recodeOrdinal(hcog, varlist_orig = varlist_orig, varlist_tr = varlist_tr)

# table(hcog$naming_hcapr)
# table(hcog$ceradtot_hcapr)
# table(hcog$ceraddel_hcapr)
# table(hcog$ceradrcg_hcapr)
# table(hcog$animfl_hcapr)
# table(hcog$bmim_hcapr)
# table(hcog$bmdel_hcapr)
# table(hcog$lmim_hcapr)
# table(hcog$lmdel_hcapr)
# table(hcog$lmrcg_hcapr)
# table(hcog$conprxim_hcapr)
# table(hcog$conprxdel_hcapr)
# table(hcog$dsmt_hcapr)
# table(hcog$numser_hcapr)
# table(hcog$ravens_hcapr)
# table(hcog$trailsa_hcapir)
# table(hcog$trailsb_hcapir)
# table(hcog$spatial_hcapr)


### End recode HCAP

saveRDS(hcog, file = "~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")
write.csv(hcog, file = "~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.csv", row.names=FALSE)


#   --------------------------------- End Data --------------------------------

Multidimensional Item Response Theory Models

#   ******************************* mirt models ********************************

hcog <- readRDS("~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")

### Item names for analysis

mmse_nmsr <- c("ornttime_mmser","orntplace_mmser","imrec_mmser","delrec_mmser",
               "spellbck_mmser","naming_mmser","phrase_mmser","commfoll_mmser","writesent_mmser",
               "draw_mmser","comm3step_mmser")

tics_nmsr <- c("wlimrc_ticsr","wldlrc_ticsr","cntbck_ticsr","animfl_ticsr",
               "numser_ticsr","serial7_ticsr","ornttime_ticsr","naming_ticsr")

hcap_nmsr <- c("naming_hcapr","ceradtot_hcapr","ceraddel_hcapr","ceradrcg_hcapr",
               "animfl_hcapr","bmim_hcapr","bmdel_hcapr","lmim_hcapr","lmdel_hcapr","lmrcg_hcapr",
               "conprxim_hcapr","conprxdel_hcapr","dsmt_hcapr","numser_hcapr","ravens_hcapr",
               "trailsa_hcapir","trailsb_hcapir","spatial_hcapr")

hcog1 <- hcog[,c(mmse_nmsr,tics_nmsr,hcap_nmsr)]
names(hcog1)

### Unidimensional model

model_hrs_hcap_mmse <- "cog = 1-37"

res_1f <- mirt(hcog1, model_hrs_hcap_mmse, 
    method="MHRM",
    technical=list(NCYCLES=1000))

summary(res_1f)

fit_1f <- M2(res_1f,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE)

### Multidimensional model - shared methods

model_hrs_hcap_mmse_m1 <- "cog = 1-37
                        recmmse = 3-4
                        wltics = 12-13
                        cerhcap = 21-23
                        bmhcap = 25-26
                        lmhcap = 27-29
                        cnpxhcap = 30-31
                        trhcap = 35-36"

res_md_1 <- mirt(hcog1, model_hrs_hcap_mmse_m1, 
               method="MHRM",
               technical=list(NCYCLES=1000))

summary(res_md_1)
coef(res_md_1)

# par <- mod2values(res_md_1)

fit_md_1 <- M2(res_md_1,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
    QMC = TRUE)

# constrained loadings of 2-indicator factors
res_md_1a <- mirt(hcog1, model_hrs_hcap_mmse_m1, 
    method="MHRM",
    constrain = list(c(28,38),c(128,145),c(310,327),c(396,412),c(479,495)),
    technical=list(NCYCLES=1000))

summary(res_md_1a)
coef(res_md_1a)

fit_md_1a <- M2(res_md_1a,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
    QMC = TRUE)


### Multidimensional model - shared methods and content

model_hrs_hcap_mmse_m2 <- "cog = 1-37
                        recmmse = 3-4
                        wltics = 12-13
                        cerhcap = 21-23
                        bmhcap = 25-26
                        lmhcap = 27-29
                        cnpxhcap = 30-31
                        trhcap = 35-36
                        ornt = 1,18
                        nmng = 6,19,20"

res_md_2 <- mirt(hcog1, model_hrs_hcap_mmse_m2, 
    method="MHRM",
    technical=list(NCYCLES=1000))

summary(res_md_2)
coef(res_md_2)

fit_md_2 <- M2(res_md_2,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
               QMC = TRUE)

# constrained loadings of 2-indicator factors
res_md_2a <- mirt(hcog1, model_hrs_hcap_mmse_m2, 
    method="MHRM",
    constrain = list(c(32,44),c(150,169),c(358,377),c(454,472),c(547,565),c(9,254)),
    technical=list(NCYCLES=1000))

# par <- mod2values(res_md_2)

summary(res_md_2a)
coef(res_md_2a)

fit_md_2a <- M2(res_md_2a,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
               QMC = TRUE)

save(res_1f,res_md_1,res_md_1a,res_md_2,res_md_2a,fit_1f,fit_md_1,fit_md_1a,fit_md_2,fit_md_2a,
     file = "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/co_calibration_results.Rdata")


# rm(list = ls()[grepl("^fit_",ls())])

#   ------------------------------end mirt models ------------------------------

Model Fit Comparisons

load("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/co_calibration_results.Rdata")

fit_1f <- fit_1f %>% mutate(model = "unidimensional") %>% 
    relocate(model)
fit_md_1a <- fit_md_1a %>% mutate(model = "multidimensional 1") %>% 
    relocate(model)
fit_md_2a <- fit_md_2a %>% mutate(model = "multidimensional 2") %>% 
    relocate(model)

fit_summ <- bind_rows(fit_1f, fit_md_1a, fit_md_2a)

pandoc.table(fit_summ, row.names = FALSE, caption = "Model fit comparisons")
Model fit comparisons (continued below)
model M2 df p RMSEA RMSEA_5 RMSEA_95 SRMSR
unidimensional 4100 629 0 0.06927 0.06723 0.07127 0.1571
multidimensional 1 1993 618 0 0.04399 0.04182 0.04614 0.1493
multidimensional 2 1991 614 0 0.04416 0.04199 0.04632 0.1488
TLI CFI
0.8226 0.8325
0.9285 0.9336
0.9279 0.9335

Test Information Curves

extractMIRTParm <- function(model,max_thresh = 9) {
    require(mirt)
    
    parm <- mod2values(model)
    
    parmw <- parm %>% dplyr::select(item, name, value) %>% 
        pivot_wider(id_cols = item, names_from = name, values_from = value) %>% 
        filter(!item == "GROUP") %>% mutate(
            d1 = case_when(
                is.na(d1) ~ d,
                TRUE ~ d1
            )
        )
    
    a <- as.matrix(parmw %>% dplyr::select(a1))
    d <- as.matrix(parmw %>% dplyr::select(paste0("d",1:max_thresh)))
    return(list(a,d))
}

hcog <- readRDS("~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")

### Item names for analysis

mmse_nmsr <- c("ornttime_mmser","orntplace_mmser","imrec_mmser","delrec_mmser",
               "spellbck_mmser","naming_mmser","phrase_mmser","commfoll_mmser","writesent_mmser",
               "draw_mmser","comm3step_mmser")

tics_nmsr <- c("wlimrc_ticsr","wldlrc_ticsr","cntbck_ticsr","animfl_ticsr",
               "numser_ticsr","serial7_ticsr","ornttime_ticsr","naming_ticsr")

hcap_nmsr <- c("naming_hcapr","ceradtot_hcapr","ceraddel_hcapr","ceradrcg_hcapr",
               "animfl_hcapr","bmim_hcapr","bmdel_hcapr","lmim_hcapr","lmdel_hcapr","lmrcg_hcapr",
               "conprxim_hcapr","conprxdel_hcapr","dsmt_hcapr","numser_hcapr","ravens_hcapr",
               "trailsa_hcapir","trailsb_hcapir","spatial_hcapr")

hcog1 <- hcog[,c(mmse_nmsr,tics_nmsr,hcap_nmsr)]

mirt_mmse <- "mmse = 1-11"
mirt_tics3 <- "tics = 1-6"
mirt_tics10 <- "tics = 1-7"
mirt_tics13 <- "tics = 1-8"
mirt_tics16 <- "tics = 1-10"
mirt_hcap <- "hcap = 1-18"
mirt_all <- "cog = 1-37"



true_sim <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/simulated_longitudinal_true_cognition.rds")

item_labels <- names(hcog1)

a <- extractMIRTParm((res_md_1a))[[1]]
d <- extractMIRTParm((res_md_1a))[[2]]

df <- data.frame(simdata(a = a, d = d, N = 1000, itemtype = 'graded'))


# cat("Test Information Curve - MMSE")
mod <- mirt(data = df[,1:11], model = mirt_mmse)
info_mmse <- plot(mod, type = 'info', ylim = c(0,40), 
     main = "Test Information Curve - MMSE")

# cat("Test Information Curve - TICS3")
mod <- mirt(data = df[,c(12:14,17:19)], model = mirt_tics3)
info_tics3 <- plot(mod, type = 'info', ylim = c(0,40), 
     main = "Test Information Curve - TICS3")

# cat("Test Information Curve - TICS10")
mod <- mirt(data = df[,c(12:14,16:19)], model = mirt_tics10)
info_tics10 <- plot(mod, type = 'info', ylim = c(0,40), 
     main = "Test Information Curve - TICS10")

# cat("Test Information Curve - TICS13")
mod <- mirt(data = df[,12:19], model = mirt_tics13)
info_tics13 <- plot(mod, type = 'info', ylim = c(0,40), 
     main = "Test Information Curve - TICS13")

# cat("Test Information Curve - TICS16")
mod <- mirt(data = df[,c(12:19,35:36)], model = mirt_tics16)
info_tics16 <- plot(mod, type = 'info', ylim = c(0,40), 
     main = "Test Information Curve - TICS16")


# cat("Test Information Curve - HCAP")
mod <- mirt(data = df[,20:37], model = mirt_hcap)
info_hcap <- plot(mod, type = 'info', ylim = c(0,40), 
     main = "Test Information Curve - HCAP")

# cat("Test Information Curve - MMSE + TICS + HCAP")
mod <- mirt(data = df, model = mirt_all)
info_all <- plot(mod, type = 'info', ylim = c(0,40), 
     main = "Test Information Curve - MMSE + TICS + HCAP")

Simulation of Longitudinal Cognitive Change

Simulation Functions

require(lme4)
require(mirt)
require(tidyverse)
require(ggplot2)

# show size of environment objects
# sort( sapply(ls(),function(x){object.size(get(x))})) 

# rm(list = ls()[grepl("^re",ls())])



#   extractMIRTParm is a function that extracts item parameters from a mirt
#       estimate multidimensional (uncorrelated dimensions) model and converts
#       discrimination (a) parameters to a vector of item parameters and 
#       threshhold (d) parameters to a matrix of d values.
#   Parameters
#       model - estimated mirt model object
#       max_thresh - maximum number of threshold paramters across items
#   Value list with 2 elements
#       a - vector of a parameters, 1 for each item
#       d - matrix of d parameters, rows correspons to items, columns to threshholds

extractMIRTParm <- function(model,max_thresh = 9) {
    require(mirt)
    
    parm <- mod2values(model)
    
    parmw <- parm %>% dplyr::select(item, name, value) %>% 
        pivot_wider(id_cols = item, names_from = name, values_from = value) %>% 
        filter(!item == "GROUP") %>% mutate(
            d1 = case_when(
                is.na(d1) ~ d,
                TRUE ~ d1
            )
        )
    
    a <- as.matrix(parmw %>% dplyr::select(a1))
    d <- as.matrix(parmw %>% dplyr::select(paste0("d",1:max_thresh)))
    return(list(a,d))
}


#   simTraj is a function that 1) estimates a mirt model using a simulated
#       item response dataset, 2) calculates factor scores based on the estimated 
#       model, 3) estimates linear mixed effects models with true cognition and  
#       simulated cognition (factor scores from simulated item responses) as dependent 
#       variables, time in study as a fixed effect variable, and person and 
#       person-by-time random effects, and 4) returns person specific estimates 
#       (random effects) of cognitive components slopes and intercepts.
#   Parameters
#       data - data frame that contains true cognition value, simulated item responses, 
#           id variable, and time in study variable
#       sim_cog_var - label for variable within dataframe (data) for which 
#           trajectories are being simulated
#       frml - formula specification for (lmer) longitudinal mixed effects model
#       iteration = iteration number passed from simulateTrajectories
#   Value - list with 3 elements:
#       data - analytic dataframe (data) with longitudinal model predicted cognition 
#           values for each assessment
#       re - data frame with estimated random effects from longitudinal model. 
#           Includes intercept and slope random effects for true cognition (_true) 
#           and simulated cognition factor scores (_sim)
#       fe - data frame with estimated fixed effects from longitudinal model. 

# simTraj <- function(data = df, item_labels = item_labels, mirt_mod = mirt_all,
#     iteration = iteration) {
simTraj <- function(data = df, sim_cog_var = "sim_cog_all", frml = frml,
    iteration = iteration) {
    require(tidyverse)
    require(lme4)
    
    # # estimate mirt model using simulated responses and generate factor scores
    re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA,
        int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA,
        int_sim_se = NA, slope_sim_se = NA)
    # 
    # tryCatch({
    #     mod <- mirt(data[,item_labels], mirt_mod)
    #     data$sim_cog <- fscores(mod)
    #     data <- data %>% relocate(sim_cog, .after = true_cog)
    # }, error=function(e){return(re_empty)})
    # data$itertion <- iteration
    
    # frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
    # time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
    # frml <- "true_cog ~ time + agebl_75 + 
    # time:agebl_75 + (1 | id)"
    # 
    tryCatch({
        # mixed effects longitudinal model - true cognition
        # res_long_true <- lmer(true_cog ~ time + (1 + time | id), data = data)
        res_long_true <- lmer(formula = frml, data = data)
        # summary(res_long_true)
        # str(coef(res_long_true))
        data$cog_true_pred <- predict(res_long_true, newdata = data, type = "response")
        re_true <- data.frame(ranef(res_long_true))
        re_true <- re_true %>% 
            pivot_wider(id_cols = grp, names_from = term,values_from = c(condval,condsd))
        names(re_true) <- c("id","int_true","slope_true","int_true_se","slope_true_se")
        re_true$id <- as.integer(levels(re_true$id))[re_true$id]
        fe_true <- data.frame(t(lme4::fixef(res_long_true)))
        names(fe_true) <- paste0(names(fe_true),"_true")
        
        
        data$sim_cog <- data[,sim_cog_var]
        #res_long_sim <- lmer(sim_cog ~ time + (1 + time | id), data = data)
        frml <- sub("true_cog","sim_cog",frml)
        res_long_sim <- lmer(formula = frml, data = data)
        # summary(res_long_sim)
        
        data$cog_sim_pred <- predict(res_long_sim, newdata = data, type = "response")
        data <- data %>% dplyr::select(-sim_cog)
        re_sim <- data.frame(ranef(res_long_sim))
        re_sim <- re_sim %>% 
            pivot_wider(id_cols = grp, names_from = term,values_from = c(condval,condsd))
        names(re_sim) <- c("id","int_sim","slope_sim","int_sim_se","slope_sim_se")
        re_sim$id <- as.integer(levels(re_sim$id))[re_sim$id] 
        fe_sim <- data.frame(t(lme4::fixef(res_long_sim)))
        names(fe_sim) <- paste0(names(fe_sim),"_sim")
        
        re <- re_true %>% left_join(re_sim, by="id")
        fe <- bind_cols(fe_true,fe_sim)
        fe$iteration <- iteration
        
    }, error=function(e){re <- re_empty})
    re$iteration <- iteration
    
    return(list(data,re,fe))
    
}

#   simulateTrajectories is a function that 1) simulates item responses for 
#   true cognition vectors for different components of the HRS-HCAP cognitive test
#   battery (MMSE, HRS TICS, HCAP, MMSE+TICS+HCAP), and 2) generates simulated observed (factor)
#   scores for each of the cognitive components (simTraj()), 3) estimates linear mixed effects models
#   with true cognition and simulated cognition (factor scores from simulated item responses)
#   as dependent variables, time in study as a fixed effect variable, and person and 
#   person-by-time random effects (simTraj()), 4) collates person specific estimates (random effects)
#   of cognitive components slopes and intercepts.
# 
#   Parameters
#       true_sim - data frame with true cognition values, rows correspond to assessments,
#           columns correspond to simulation iterations.
#       a_par - vector/data frame containing a (discrimination) item parameters 
#           from mirt co-calibration model
#       d_par - vector/data frame containing d (threshhold) item parameters 
#           from mirt co-calibration model
#       item_labels - item labels
#       iter_group - number of groups of (niter) iterations
#       niter - number of simulation iterations within iteration groups
#       out_dir - path to directory for output of simulation results
#       frml - formula specification for (lmer) longitudinal mixed effects model
#
#   Value - Saves lists of results to out_dir - there is a file for each iter_group 
#       that contains niter elements:
#       ds - list of datasets for each niter simulations. Each dataset contains
#           the true cognition value (true_cog), the simulated item responses conditional
#           on true_cog, simulated item response data and cognition factor scores, 
#           and longitudinal mixed effect model predicted cognition values for each simulated
#           assessment. Measures include MMSE, TICS, HCAP, and MMSE+TICS+HCAP.
#       re_all - estimated random effects for the full HRS-HCAP cognitive test battery
#           (MMSE+TICS+HCAP). Includes intercept and slope random effects for true
#           cognition (_true) and simulated cognition factor scores (_sim)
#       re_mmse - estimated random effects for the MMSE component of the HRS-HCAP 
#           cognitive test battery. Includes intercept and slope random effects for true
#           cognition (_true) and simulated cognition factor scores (_sim)
#       re_tics - estimated random effects for the TICS component of the HRS-HCAP 
#           cognitive test battery. Includes intercept and slope random effects for true
#           cognition (_true) and simulated cognition factor scores (_sim)
#       re_hcap - estimated random effects for the HCAP component of the HRS-HCAP 
#           cognitive test battery. Includes intercept and slope random effects for true
#           cognition (_true) and simulated cognition factor scores (_sim)
#       fe_all - estimated fixed effects for mixed effect model of cognition measured 
#           by full HRS-HCAP cognitive test battery (MMSE+TICS+HCAP). 
#           Includes intercept and slope random effects for true
#           cognition (_true) and simulated cognition factor scores (_sim).
#       fe_mmse - estimated fixed effects for mixed effect model of cognition measured 
#           by MMSE. 
#           Includes intercept and slope random effects for true
#           cognition (_true) and simulated cognition factor scores (_sim).
#       fe_tics - estimated fixed effects for mixed effect model of cognition measured 
#           by TICS. 
#           Includes intercept and slope random effects for true
#           cognition (_true) and simulated cognition factor scores (_sim).
#       fe_hcap - estimated fixed effects for mixed effect model of cognition measured 
#           by HCAP. 
#           Includes intercept and slope random effects for true
#           cognition (_true) and simulated cognition factor scores (_sim).

simulateTrajectories <- function(true_sim, a_par, d_par, item_labels, 
    niter = 20, iter_group = 5, out_dir = "Analysis/Simulation Results/",
    frml = "true_cog ~ time + (1 | id)") {
    require(mirt)
    require(tidyverse)
    for (itgrp in 1:iter_group) {
        ds <- list()
        re_all <- list()
        re_mmse <- list()
        re_tics <- list()
        re_hcap <- list()
        fe_all <- list()
        fe_mmse <- list()
        fe_tics <- list()
        fe_hcap <- list()
        set.seed(092724)
        for (iter in 1:niter) {
            cat(paste0("Group - ",itgrp, ", Iteration - ",iter,"\n"))
            
            # simulate item responses
            iteration <- ((itgrp - 1) * niter) + iter
            df <- data.frame(simdata(a = a_par, d = d_par, Theta = true_sim[,paste0("vm",iteration)], 
                itemtype = 'graded'))
            names(df) <- item_labels
            df$id <- true_sim$id
            df$time <- true_sim$time
            df$agebl_75 <- true_sim$agebl_75
            df$true_cog <- true_sim[,paste0("vm",iteration)]
            df$iteration <- iteration
            df <- df %>% relocate(c(id,iteration,time,true_cog))
            
            # true random effect intercept and slope - from calculate_re_true.R 
            true_re <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/true_random_effects.rds")
            df <- df %>% left_join((true_re %>% dplyr::select(id,int_true_qrtl,
                slope_true_qrtl)),by="id")
            
            flag_low_response <- FALSE
            for (itnm in item_labels) {
                if (length(table(df[,itnm])) < 2) {
                    flag_low_response <- TRUE
                }
            }
            if (flag_low_response == TRUE) {
                iter = iter - 1
                next
            }
            
            # mirt_models for cognition measures
            
            mirt_mmse <- "mmse = 1-11"
            mirt_tics <- "tics = 1-8"
            mirt_hcap <- "hcap = 1-18"
            mirt_all <- "cog = 1-37"
            
            ### estimate mirt model using simulated responses and generate factor scores
            # create empty re dataframe to return if error in generation
            re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA, 
                int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA, 
                int_sim_se = NA, slope_sim_se = NA)
            
            tryCatch({
                mod <- mirt(df[,item_labels], mirt_all)
                df$sim_cog_all <- fscores(mod)
                mod <- mirt(df[,item_labels[1:11]], mirt_mmse)
                df$sim_cog_mmse <- fscores(mod)
                mod <- mirt(df[,item_labels[12:19]], mirt_tics)
                df$sim_cog_tics <- fscores(mod)
                mod <- mirt(df[,item_labels[20:37]], mirt_hcap)
                df$sim_cog_hcap <- fscores(mod)
                df <- df %>% 
                    relocate(c(sim_cog_all,sim_cog_mmse,sim_cog_tics,sim_cog_hcap), 
                        .after = true_cog)
            }, error=function(e){return(re_empty)})
            
            # rescale factor scores to equate metric to true cognition
            df <- df %>% mutate(
                sim_cog_mmse_rs = (((sim_cog_mmse - mean(sim_cog_mmse)) / 
                    sd(sim_cog_mmse)) * sd(true_cog)) + mean(true_cog),
                sim_cog_tics_rs = (((sim_cog_tics - mean(sim_cog_tics)) / 
                    sd(sim_cog_tics)) * sd(true_cog)) + mean(true_cog),
                sim_cog_hcap_rs = (((sim_cog_hcap - mean(sim_cog_hcap)) / 
                    sd(sim_cog_hcap)) * sd(true_cog)) + mean(true_cog),
                sim_cog_all_rs = (((sim_cog_all - mean(sim_cog_all)) / 
                    sd(sim_cog_all)) * sd(true_cog)) + mean(true_cog),
            )
            
            ### calculate simulated random effects
            
            # frml <- as.formula(frml)
            # frml <- as.formula("true_cog ~ time + (1 | id)")
            # frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
            #     time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
            
            # MMSE+TICS+HCAP
            res <- simTraj(data = df, sim_cog_var = "sim_cog_all_rs",frml = frml,
                iteration = iteration)
            df <- res[[1]]
            df <- df %>% rename(cog_true_pred_all = cog_true_pred,
                cog_sim_pred_all = cog_sim_pred)
            re <- res[[2]]
            re$label <- "MMSE+TICS+HCAP"
            re_all[[paste0("iteration-",iteration)]] <- re
            fe <- res[[3]]
            fe$label <- "MMSE+TICS+HCAP"
            fe_all[[paste0("iteration-",iteration)]] <- fe
            
            # MMSE
            res <- simTraj(data = df, sim_cog_var = "sim_cog_mmse_rs",frml = frml,
                iteration = iteration)
            df <- res[[1]]
            df <- df %>% rename(cog_true_pred_mmse = cog_true_pred,
                 cog_sim_pred_mmse = cog_sim_pred)
            re <- res[[2]]
            re$label <- "MMSE"
            re_mmse[[paste0("iteration-",iteration)]] <- re
            fe <- res[[3]]
            fe$label <- "MMSE"
            fe_mmse[[paste0("iteration-",iteration)]] <- fe
            
            # TICS
            res <- simTraj(data = df, sim_cog_var = "sim_cog_tics_rs",frml = frml,
                iteration = iteration)
            df <- res[[1]]
            df <- df %>% rename(cog_true_pred_tics = cog_true_pred,
                cog_sim_pred_tics = cog_sim_pred)
            re <- res[[2]]
            re$label <- "TICS"
            re_tics[[paste0("iteration-",iteration)]] <- re
            fe <- res[[3]]
            fe$label <- "TICS"
            fe_tics[[paste0("iteration-",iteration)]] <- fe
            
            # HCAP
            res <- simTraj(data = df, sim_cog_var = "sim_cog_hcap_rs",frml = frml,
                iteration = iteration)
            df <- res[[1]]
            df <- df %>% rename(cog_true_pred_hcap = cog_true_pred,
                cog_sim_pred_hcap = cog_sim_pred)
            re <- res[[2]]
            re$label <- "HCAP"
            re_hcap[[paste0("iteration-",iteration)]] <- re
            fe <- res[[3]]
            fe$label <- "HCAP"
            fe_hcap[[paste0("iteration-",iteration)]] <- fe
            
            ds[[paste0("iteration-",iteration)]] <- df
            
            
        } # end for iter
        
        saveRDS(ds, file = paste0(out_dir,"ds_iteration_group_", itgrp, ".rds"))
        saveRDS(re_mmse, file = paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
        saveRDS(re_tics, file = paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
        saveRDS(re_hcap, file = paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
        saveRDS(re_all, file = paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
        saveRDS(fe_mmse, file = paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
        saveRDS(fe_tics, file = paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
        saveRDS(fe_hcap, file = paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
        saveRDS(fe_all, file = paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
        
    } # end for itgrp
    
    for (itgrp in 1:iter_group) {
        if (itgrp == 1) {
            re_mmse <- readRDS(paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
            re_mmse <- re_mmse %>% bind_rows()
            re_tics <- readRDS(paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
            re_tics <- re_tics %>% bind_rows()
            re_hcap <- readRDS(paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
            re_hcap <- re_hcap %>% bind_rows()
            re_all <- readRDS(paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
            re_all <- re_all %>% bind_rows()
            
            fe_mmse <- readRDS(paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
            fe_mmse <- fe_mmse %>% bind_rows()
            fe_tics <- readRDS(paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
            fe_tics <- fe_tics %>% bind_rows()
            fe_hcap <- readRDS(paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
            fe_hcap <- fe_hcap %>% bind_rows()
            fe_all <- readRDS(paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
            fe_all <- fe_all %>% bind_rows()
            
        } else {
            re_mmse_temp <- readRDS(paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
            re_mmse_temp <- re_mmse_temp %>% bind_rows()
            re_mmse <- bind_rows(re_mmse, re_mmse_temp)
            re_tics_temp <- readRDS(paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
            re_tics_temp <- re_tics_temp %>% bind_rows()
            re_tics <- bind_rows(re_tics, re_tics_temp)
            re_hcap_temp <- readRDS(paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
            re_hcap_temp <- re_hcap_temp %>% bind_rows()
            re_hcap <- bind_rows(re_hcap, re_hcap_temp)
            re_all_temp <- readRDS(paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
            re_all_temp <- re_all_temp %>% bind_rows()
            re_all <- bind_rows(re_all, re_all_temp)
            
            fe_mmse_temp <- readRDS(paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
            fe_mmse_temp <- fe_mmse_temp %>% bind_rows()
            fe_mmse <- bind_rows(fe_mmse, fe_mmse_temp)
            fe_tics_temp <- readRDS(paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
            fe_tics_temp <- fe_tics_temp %>% bind_rows()
            fe_tics <- bind_rows(fe_tics, fe_tics_temp)
            fe_hcap_temp <- readRDS(paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
            fe_hcap_temp <- fe_hcap_temp %>% bind_rows()
            fe_hcap <- bind_rows(fe_hcap, fe_hcap_temp)
            fe_all_temp <- readRDS(paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
            fe_all_temp <- fe_all_temp %>% bind_rows()
            fe_all <- bind_rows(fe_all, fe_all_temp)
            
        }
    }
    saveRDS(re_mmse, file = paste0(out_dir,"re_mmse", ".rds"))
    saveRDS(re_tics, file = paste0(out_dir,"re_tics", ".rds"))
    saveRDS(re_hcap, file = paste0(out_dir,"re_hcap", ".rds"))
    saveRDS(re_all, file = paste0(out_dir,"re_all", ".rds"))

    saveRDS(fe_mmse, file = paste0(out_dir,"fe_mmse", ".rds"))
    saveRDS(fe_tics, file = paste0(out_dir,"fe_tics", ".rds"))
    saveRDS(fe_hcap, file = paste0(out_dir,"fe_hcap", ".rds"))
    saveRDS(fe_all, file = paste0(out_dir,"fe_all", ".rds"))
    
    
    # return(list("ds" = ds, "re_all" = re_all, "re_mmse" = re_mmse,
    #     "re_tics" = re_tics, "re_hcap" = re_hcap))
}


simulateTrajectories2 <- function(item_labels, niter = 20, iter_group = 5,
        in_dir = "Analysis/Simulation Results/",
        out_dir = "Analysis/Simulation Results/",
        frml = "true_cog ~ time + (1 | id)") {
    require(mirt)
    require(tidyverse)
    for (itgrp in 1:iter_group) {
        ds <- list()
        re_tics3 <- list()
        re_tics10 <- list()
        re_tics13 <- list()
        re_tics16 <- list()
        fe_tics3 <- list()
        fe_tics10 <- list()
        fe_tics13 <- list()
        fe_tics16 <- list()
        set.seed(092724)
        for (iter in 1:niter) {
            cat(paste0("Group - ",itgrp, ", Iteration - ",iter,"\n"))
            
            #  load simulated item responses
            iteration <- ((itgrp - 1) * niter) + iter
            df <- readRDS(paste0(in_dir,"ds_iteration_group_",itgrp,".rds"))[[iter]]
            df <- df %>% dplyr::select(c(id:true_cog,agebl_75:slope_true_qrtl,any_of(item_labels)))
            
            # ds_iteration_group_36.rds
            # 
            # df <- data.frame(simdata(a = a_par, d = d_par, Theta = true_sim[,paste0("vm",iteration)], 
            #                          itemtype = 'graded'))
            # names(df) <- item_labels
            # df$id <- true_sim$id
            # df$time <- true_sim$time
            # df$agebl_75 <- true_sim$agebl_75
            # df$true_cog <- true_sim[,paste0("vm",iteration)]
            # df$iteration <- iteration
            # df <- df %>% relocate(c(id,iteration,time,true_cog))
            
            # # true random effect intercept and slope - from calculate_re_true.R 
            # true_re <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/true_random_effects.rds")
            # df <- df %>% left_join((true_re %>% dplyr::select(id,int_true_qrtl,
            #                                                   slope_true_qrtl)),by="id")
            # 
            # flag_low_response <- FALSE
            # for (itnm in item_labels) {
            #     if (length(table(df[,itnm])) < 2) {
            #         flag_low_response <- TRUE
            #     }
            # }
            # if (flag_low_response == TRUE) {
            #     iter = iter - 1
            #     next
            # }
            
            # mirt_models for cognition measures
            
            # mirt_mmse <- "mmse = 1-11"
            # mirt_tics <- "tics = 1-8"
            # mirt_hcap <- "hcap = 1-18"
            # mirt_all <- "cog = 1-37"

            mirt_tics3 <- "tics3 = 1-6"
            mirt_tics10 <- "tics10 = 1-7"
            mirt_tics13 <- "tics13 = 1-8"
            mirt_tics16 <- "tics16 = 1-10"
            
            ### estimate mirt model using simulated responses and generate factor scores
            # create empty re dataframe to return if error in generation
            re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA, 
                                   int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA, 
                                   int_sim_se = NA, slope_sim_se = NA)
            
            tryCatch({
                mod <- mirt(df[,item_labels[c(1:3,6:8)]], mirt_tics3)
                df$sim_cog_tics3 <- fscores(mod)
                mod <- mirt(df[,item_labels[c(1:3,5:8)]], mirt_tics10)
                df$sim_cog_tics10 <- fscores(mod)
                mod <- mirt(df[,item_labels[1:8]], mirt_tics13)
                df$sim_cog_tics13 <- fscores(mod)
                mod <- mirt(df[,item_labels[1:10]], mirt_tics16)
                df$sim_cog_tics16 <- fscores(mod)
                df <- df %>% 
                    relocate(c(sim_cog_tics3,sim_cog_tics10,sim_cog_tics13,sim_cog_tics16), 
                        .after = true_cog)
            }, error=function(e){return(re_empty)})
            
            # rescale factor scores to equate metric to true cognition
            df <- df %>% mutate(
                sim_cog_tics3_rs = (((sim_cog_tics3 - mean(sim_cog_tics3)) / 
                    sd(sim_cog_tics3)) * sd(true_cog)) + mean(true_cog),
                sim_cog_tics10_rs = (((sim_cog_tics10 - mean(sim_cog_tics10)) / 
                    sd(sim_cog_tics10)) * sd(true_cog)) + mean(true_cog),
                sim_cog_tics13_rs = (((sim_cog_tics13 - mean(sim_cog_tics13)) / 
                    sd(sim_cog_tics13)) * sd(true_cog)) + mean(true_cog),
                sim_cog_tics16_rs = (((sim_cog_tics16 - mean(sim_cog_tics16)) / 
                    sd(sim_cog_tics16)) * sd(true_cog)) + mean(true_cog),
            )
            
            ### calculate simulated random effects
            
            # frml <- as.formula(frml)
            # frml <- as.formula("true_cog ~ time + (1 | id)")
            # frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
            #     time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
            
            # TICS3
            res <- simTraj(data = df, sim_cog_var = "sim_cog_tics3_rs",frml = frml,
                iteration = iteration)
            df <- res[[1]]
            df <- df %>% rename(cog_true_pred_tics3 = cog_true_pred,
                cog_sim_pred_tics3 = cog_sim_pred)
            re <- res[[2]]
            re$label <- "TICS3"
            re_tics3[[paste0("iteration-",iteration)]] <- re
            fe <- res[[3]]
            fe$label <- "TICS3"
            fe_tics3[[paste0("iteration-",iteration)]] <- fe
            
            # TICS10
            res <- simTraj(data = df, sim_cog_var = "sim_cog_tics10_rs",frml = frml,
                iteration = iteration)
            df <- res[[1]]
            df <- df %>% rename(cog_true_pred_tics10 = cog_true_pred,
                cog_sim_pred_tics10 = cog_sim_pred)
            re <- res[[2]]
            re$label <- "TICS10"
            re_tics10[[paste0("iteration-",iteration)]] <- re
            fe <- res[[3]]
            fe$label <- "TICS10"
            fe_tics10[[paste0("iteration-",iteration)]] <- fe
            
            # TICS13
            res <- simTraj(data = df, sim_cog_var = "sim_cog_tics13_rs",frml = frml,
                           iteration = iteration)
            df <- res[[1]]
            df <- df %>% rename(cog_true_pred_tics13 = cog_true_pred,
                cog_sim_pred_tics13 = cog_sim_pred)
            re <- res[[2]]
            re$label <- "TICS13"
            re_tics13[[paste0("iteration-",iteration)]] <- re
            fe <- res[[3]]
            fe$label <- "TICS13"
            fe_tics13[[paste0("iteration-",iteration)]] <- fe
            
            # TICS16
            res <- simTraj(data = df, sim_cog_var = "sim_cog_tics16_rs",frml = frml,
                iteration = iteration)
            df <- res[[1]]
            df <- df %>% rename(cog_true_pred_tics16 = cog_true_pred,
                cog_sim_pred_tics16 = cog_sim_pred)
            re <- res[[2]]
            re$label <- "TICS16"
            re_tics16[[paste0("iteration-",iteration)]] <- re
            fe <- res[[3]]
            fe$label <- "TICS16"
            fe_tics16[[paste0("iteration-",iteration)]] <- fe
            
            ds[[paste0("iteration-",iteration)]] <- df
            
            
        } # end for iter
        
        saveRDS(ds, file = paste0(out_dir,"ds_iteration_group_", itgrp, ".rds"))
        saveRDS(re_tics3, file = paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
        saveRDS(re_tics10, file = paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
        saveRDS(re_tics13, file = paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
        saveRDS(re_tics16, file = paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
        saveRDS(fe_tics3, file = paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
        saveRDS(fe_tics10, file = paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
        saveRDS(fe_tics13, file = paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
        saveRDS(fe_tics16, file = paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
        
    } # end for itgrp
    
    for (itgrp in 1:iter_group) {
        if (itgrp == 1) {
            re_tics3 <- readRDS(paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
            re_tics3 <- re_tics3 %>% bind_rows()
            re_tics10 <- readRDS(paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
            re_tics10 <- re_tics10 %>% bind_rows()
            re_tics13 <- readRDS(paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
            re_tics13 <- re_tics13 %>% bind_rows()
            re_tics16 <- readRDS(paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
            re_tics16 <- re_tics16 %>% bind_rows()
            
            fe_tics3 <- readRDS(paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
            fe_tics3 <- fe_tics3 %>% bind_rows()
            fe_tics10 <- readRDS(paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
            fe_tics10 <- fe_tics10 %>% bind_rows()
            fe_tics13 <- readRDS(paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
            fe_tics13 <- fe_tics13 %>% bind_rows()
            fe_tics16 <- readRDS(paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
            fe_tics16 <- fe_tics16 %>% bind_rows()
            
        } else {
            re_tics3_temp <- readRDS(paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
            re_tics3_temp <- re_tics3_temp %>% bind_rows()
            re_tics3 <- bind_rows(re_tics3, re_tics3_temp)
            re_tics10_temp <- readRDS(paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
            re_tics10_temp <- re_tics10_temp %>% bind_rows()
            re_tics10 <- bind_rows(re_tics10, re_tics10_temp)
            re_tics13_temp <- readRDS(paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
            re_tics13_temp <- re_tics13_temp %>% bind_rows()
            re_tics13 <- bind_rows(re_tics13, re_tics13_temp)
            re_tics16_temp <- readRDS(paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
            re_tics16_temp <- re_tics16_temp %>% bind_rows()
            re_tics16 <- bind_rows(re_tics16, re_tics16_temp)
            
            fe_tics3_temp <- readRDS(paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
            fe_tics3_temp <- fe_tics3_temp %>% bind_rows()
            fe_tics3 <- bind_rows(fe_tics3, fe_tics3_temp)
            fe_tics10_temp <- readRDS(paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
            fe_tics10_temp <- fe_tics10_temp %>% bind_rows()
            fe_tics10 <- bind_rows(fe_tics10, fe_tics10_temp)
            fe_tics13_temp <- readRDS(paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
            fe_tics13_temp <- fe_tics13_temp %>% bind_rows()
            fe_tics13 <- bind_rows(fe_tics13, fe_tics13_temp)
            fe_tics16_temp <- readRDS(paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
            fe_tics16_temp <- fe_tics16_temp %>% bind_rows()
            fe_tics16 <- bind_rows(fe_tics16, fe_tics16_temp)
            
        }
    }
    saveRDS(re_tics3, file = paste0(out_dir,"re_tics3", ".rds"))
    saveRDS(re_tics10, file = paste0(out_dir,"re_tics10", ".rds"))
    saveRDS(re_tics13, file = paste0(out_dir,"re_tics13", ".rds"))
    saveRDS(re_tics16, file = paste0(out_dir,"re_tics16", ".rds"))
    
    saveRDS(fe_tics3, file = paste0(out_dir,"fe_tics3", ".rds"))
    saveRDS(fe_tics10, file = paste0(out_dir,"fe_tics10", ".rds"))
    saveRDS(fe_tics13, file = paste0(out_dir,"fe_tics13", ".rds"))
    saveRDS(fe_tics16, file = paste0(out_dir,"fe_tics16", ".rds"))
    
    
    # return(list("ds" = ds, "re_all" = re_all, "re_mmse" = re_mmse,
    #     "re_tics" = re_tics, "re_hcap" = re_hcap))
}


# dat_cols <- c("id","time","age","agebl_75","slope_true_qrtl","int_true_qrtl")
# res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
#               "cog_sim_pred_tics","cog_sim_pred_hcap")

mergeSimResults <-function(file_name = "ds_iteration_group_1.rds",
                           out_dir = out_dir,
                           dat_cols = dat_cols, 
                           res_cols = res_cols,
                           iter_group = 50,
                           niter = 20) {
    res <- (readRDS(paste0(out_dir,file_name))[[1]] %>%
                dplyr::select(any_of(dat_cols)))
    for (i in 1:iter_group) {
        for (n in 1:niter) {
            iteration <- ((i-1) * niter) + n
            fl_nm <- sub("_1",paste0("_",i),file_name)
            ds <- readRDS(paste0(out_dir,fl_nm))[[n]] %>% 
                dplyr::select(all_of(res_cols)) 
            names(ds) <- paste0(res_cols,".",iteration)
            res <- res %>% bind_cols(ds)
        }
    }   
        return(res)
}

Simulated Cognition Change versus True Cognition Change Plots

# sim20 <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/sim20.rds")
# 
# re_all <- sim20[["re_all"]]
# re_mmse <- sim20[["re_mmse"]]
# re_tics <- sim20[["re_tics"]]
# re_hcap <- sim20[["re_hcap"]]

out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/10_04_24_unconditional/"

re_all <- readRDS(paste0(out_dir,"re_all.rds"))
re_mmse <- readRDS(paste0(out_dir,"re_mmse.rds"))
re_tics <- readRDS(paste0(out_dir,"re_tics.rds"))
re_hcap <- readRDS(paste0(out_dir,"re_hcap.rds"))

re_all_summ <- bind_rows(re_all)
re_mmse_summ <- bind_rows(re_mmse)
re_tics_summ <- bind_rows(re_tics)
re_hcap_summ <- bind_rows(re_hcap)

ggplot(data=re_mmse,aes(slope_true,slope_sim)) + 
    geom_point() + 
    geom_smooth() +
    labs(caption = "True versus simulated cognitive change - MMSE")

ggplot(data=re_tics,aes(slope_true,slope_sim)) + 
    geom_point() + 
    geom_smooth() +
    labs(caption = "True versus simulated cognitive change - TICS")

ggplot(data=re_hcap,aes(slope_true,slope_sim)) + 
    geom_point() + 
    geom_smooth() +
    labs(caption = "True versus simulated cognitive change - HCAP")

ggplot(data=re_all,aes(slope_true,slope_sim)) + 
    geom_point() + 
    geom_smooth() +
    labs(caption = "True versus simulated cognitive change - MMSE + TICS + HCAP")

Simulated Cognitive Trajectories versus True Cognition Trajectories

# Create model predicted trajectories 

# # code to create merged r1 dataset
# dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
# res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
#               "cog_sim_pred_tics","cog_sim_pred_hcap")
# out_dir <- "Analysis/Simulation Results/"
# 
# r1 <- mergeSimResults(file_name = "ds_iteration_group_1.rds", out_dir = out_dir, 
#                       dat_cols = dat_cols, res_cols = res_cols, iter_group = 50, niter = 20)
# saveRDS(r1,file=paste0(out_dir,"model_predicted_trajectory_simulations.rds"))
# 
# rm(r1)


# prepare longitudinal data file

out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/10_15_24/"
dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
              "cog_sim_pred_tics","cog_sim_pred_hcap")


r1 <- readRDS(paste0(out_dir,"model_predicted_trajectory_simulations.rds"))

r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_all_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_all_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_mmse_mn <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_se <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_tics_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics_se <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_hcap_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>% 
    mutate(age = agebl_75 + 75)
# summary(r2$age)

r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_hcap_se)) %>% 
    pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_hcap_mn,
        names_to = "label", values_to = "cog") %>% mutate(
          label = sub("_mn$","",label)
        )
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_hcap_mn)) %>% 
    pivot_longer(cols = cog_true_pred_se:cog_sim_pred_hcap_se,
        names_to = "label", values_to = "cog_se") %>% mutate(
          label = sub("_se$","",label)
        )

r3 <- r3a %>% left_join((r3b %>% 
    dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>% 
    mutate(
        cognitive_measure = factor(case_when(
            label == "cog_true_pred" ~ "True Cognition",
            label == "cog_sim_pred_all" ~ "MMSE+TICS+HCAP",
            label == "cog_sim_pred_mmse" ~ "MMSE",
            label == "cog_sim_pred_tics" ~ "TICS",
            label == "cog_sim_pred_hcap" ~ "HCAP",
        ), levels = c("True Cognition","MMSE+TICS+HCAP","HCAP","TICS","MMSE"))
    )

r4 <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
    mean_cog = mean(cog),
    se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
    lower = mean_cog - (1.96 * se_cog),
    upper = mean_cog + (1.96 * se_cog),
)

r2ta <- r2 %>% dplyr::select(id,time,agebl_75,slope_true_qrtl,
    cog_sim_pred_tics_mn,cog_sim_pred_all_mn,cog_sim_pred_tics_se,cog_sim_pred_all_se) %>%
      mutate(
        cog_blend_mn = case_when(
          time <= 5 ~ cog_sim_pred_tics_mn,
          time > 5 ~ cog_sim_pred_all_mn
        ),
        cog_blend_se = case_when(
          time <= 5 ~ cog_sim_pred_tics_se,
          time > 5 ~ cog_sim_pred_all_se
        )
      ) %>% relocate(cog_blend_mn, .after = cog_sim_pred_all_mn)

r3taa <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics_se:cog_blend_se,cog_sim_pred_all_mn)) %>% 
    pivot_longer(cols = c(cog_sim_pred_tics_mn,cog_blend_mn),
        names_to = "label", values_to = "cog") %>% mutate(
          label = sub("_mn$","",label)
        )
r3tab <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics_mn:cog_blend_mn,cog_sim_pred_all_se)) %>% 
    pivot_longer(cols = c(cog_sim_pred_tics_se,cog_blend_se),
        names_to = "label", values_to = "cog_se") %>% mutate(
          label = sub("_se$","",label)
        )

r3ta <- r3taa %>% left_join((r3tab %>% 
    dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>% 
    mutate(
        cognitive_measure = factor(case_when(
            label == "cog_sim_pred_tics" ~ "TICS",
            label == "cog_blend" ~ "Blended TICS,MMSE+TICS+HCAP"
        ), levels = c("Blended TICS,MMSE+TICS+HCAP","TICS"))
    )

r4ta <- r3ta %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
    mean_cog = mean(cog),
    se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
    lower = mean_cog - (1.96 * se_cog),
    upper = mean_cog + (1.96 * se_cog),
)

rm(r2,r3,r3a,r3b,r2ta,r3ta,r3taa,r3tab)

# cor(r1[,6:31])

### Repeat for iteration 1:100

r1 <- r1[,1:505]


r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_all_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_all_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_mmse_mn <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_se <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_tics_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics_se <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_hcap_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>% 
    mutate(age = agebl_75 + 75)
# summary(r2$age)

r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_hcap_se)) %>% 
    pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_hcap_mn,
        names_to = "label", values_to = "cog") %>% mutate(
          label = sub("_mn$","",label)
        )
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_hcap_mn)) %>% 
    pivot_longer(cols = cog_true_pred_se:cog_sim_pred_hcap_se,
        names_to = "label", values_to = "cog_se") %>% mutate(
          label = sub("_se$","",label)
        )

r3 <- r3a %>% left_join((r3b %>% 
    dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>% 
    mutate(
        cognitive_measure = factor(case_when(
            label == "cog_true_pred" ~ "True Cognition",
            label == "cog_sim_pred_all" ~ "MMSE+TICS+HCAP",
            label == "cog_sim_pred_mmse" ~ "MMSE",
            label == "cog_sim_pred_tics" ~ "TICS",
            label == "cog_sim_pred_hcap" ~ "HCAP",
        ), levels = c("True Cognition","MMSE+TICS+HCAP","HCAP","TICS","MMSE"))
    )

r4a <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
    mean_cog = mean(cog),
    se_cog = mean(cog_se) / sqrt(100)
) %>% mutate(
    lower = mean_cog - (1.96 * se_cog),
    upper = mean_cog + (1.96 * se_cog),
)

rm(r1,r2,r3,r3a,r3b)




### End data   

pd <- position_dodge(width = 0.2) # move them .2 to the left and right

ggplot(data = r4, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    facet_wrap(vars(cognitive_measure)) +
    labs(caption = "True versus simulated cognitive trajectories (1000 Simulations) - All Cognition Measures")

ggplot(data = r4a, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    facet_wrap(vars(cognitive_measure)) +
    labs(caption = "True versus simulated cognitive trajectories (100 Simulations) - All Cognition Measures")

r4 <- r4 %>% mutate(n_sim = 1000)
r4a <- r4a %>% mutate(n_sim = 100)


r10 <- r4 %>% bind_rows(r4a) %>% mutate(
  n_sim = factor(n_sim, levels = c(100,1000))
)

ggplot(data = r10, aes(x = time, y = mean_cog, color = slope_true_qrtl,
        linetype = n_sim)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    facet_wrap(vars(cognitive_measure)) +
    labs(caption = "True versus simulated cognitive trajectories by number of simulations - All Cognition Measures")

r5 <- r4 %>% filter(cognitive_measure %in% c("MMSE+TICS+HCAP","True Cognition"))

ggplot(data = r5, aes(x = time, y = mean_cog, color = slope_true_qrtl,
    linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = "True versus simulated cognitive change (1000 Simulations) - MMSE + TICS + HCAP")

r6 <- r4 %>% filter(cognitive_measure %in% c("HCAP","True Cognition"))

ggplot(data = r6, aes(x = time, y = mean_cog, color = slope_true_qrtl,
    linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = "True versus simulated cognitive change (1000 Simulations) - HCAP")

r7 <- r4 %>% filter(cognitive_measure %in% c("TICS","True Cognition"))

ggplot(data = r7, aes(x = time, y = mean_cog, color = slope_true_qrtl,
    linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = "True versus simulated cognitive change (1000 Simulations) - TICS")

r7a <- r4a %>% filter(cognitive_measure %in% c("TICS","True Cognition"))

ggplot(data = r7a, aes(x = time, y = mean_cog, color = slope_true_qrtl,
    linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = "True versus simulated cognitive change (100 Simulations) - TICS")

r8 <- r4 %>% filter(cognitive_measure %in% c("MMSE","True Cognition"))

ggplot(data = r8, aes(x = time, y = mean_cog, color = slope_true_qrtl,
    linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = "True versus simulated cognitive change (1000 Simulations) - MMSE")

r8a <- r4a %>% filter(cognitive_measure %in% c("MMSE","True Cognition"))

ggplot(data = r8a, aes(x = time, y = mean_cog, color = slope_true_qrtl,
    linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = "True versus simulated cognitive change (100 Simulations) - MMSE")

r9 <- r4 %>% filter(cognitive_measure %in% c("MMSE+TICS+HCAP","MMSE"))

ggplot(data = r9, aes(x = time, y = mean_cog, color = slope_true_qrtl,
        linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = "Simulated cognitive change (1000 Simulations) - MMSE + TICS + HCAP vs. MMSE")

ggplot(data = r4ta, aes(x = time, y = mean_cog, color = slope_true_qrtl,
      linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = str_wrap("TICS all assessments versus TICS (assessments bl-5) blended with MMSE+TICS+HCAP (assessments 6-10) (1000 Simulations)",width = 70))

rm(r4,r4a,r5,r6,r7,r7a,r8,r8a,r9,r10)

Simulated Cognitive Trajectories versus True Cognition Trajectories - TICS Variants

out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/tics/"
dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
res_cols <- c("cog_true_pred_tics3","cog_sim_pred_tics3","cog_sim_pred_tics10",
              "cog_sim_pred_tics13","cog_sim_pred_tics16")


r1 <- readRDS(paste0(out_dir,"model_predicted_trajectory_simulations.rds"))

r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_tics3_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics3", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics3_se <- apply(r1[,names(r1)[grepl("sim_pred_tics3", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_tics10_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics10", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics10_se <- apply(r1[,names(r1)[grepl("sim_pred_tics10", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_tics13_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics13", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics13_se <- apply(r1[,names(r1)[grepl("sim_pred_tics13", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r1$cog_sim_pred_tics16_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics16", names(r1))]], 1,
    function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics16_se <- apply(r1[,names(r1)[grepl("sim_pred_tics16", names(r1))]], 1,
    function(x) sd(x, na.rm = TRUE))

r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>% 
    mutate(age = agebl_75 + 75)
# summary(r2$age)

r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_tics16_se)) %>% 
    pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_tics16_mn,
        names_to = "label", values_to = "cog") %>% mutate(
          label = sub("_mn$","",label)
        )
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_tics16_mn)) %>% 
    pivot_longer(cols = cog_true_pred_se:cog_sim_pred_tics16_se,
        names_to = "label", values_to = "cog_se") %>% mutate(
          label = sub("_se$","",label)
        )

r3 <- r3a %>% left_join((r3b %>% 
    dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>% 
    mutate(
        cognitive_measure = factor(case_when(
            label == "cog_true_pred" ~ "True Cognition",
            label == "cog_sim_pred_tics3" ~ "TICS (wave 3)",
            label == "cog_sim_pred_tics10" ~ "TICS (wave 3 + Number Series)",
            label == "cog_sim_pred_tics13" ~ "TICS (wave 13)",
            label == "cog_sim_pred_tics16" ~ "TICS (wave 16)"
        ), levels = c("True Cognition","TICS (wave 16)","TICS (wave 13)",
            "TICS (wave 3 + Number Series)","TICS (wave 3)"))
    )

r4 <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
    mean_cog = mean(cog),
    se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
    lower = mean_cog - (1.96 * se_cog),
    upper = mean_cog + (1.96 * se_cog),
)

r2ta <- r2 %>% dplyr::select(id,time,agebl_75,slope_true_qrtl,
    cog_sim_pred_tics3_mn,cog_sim_pred_tics16_mn,cog_sim_pred_tics3_se,cog_sim_pred_tics16_se) %>%
      mutate(
        cog_blend_mn = case_when(
          time <= 5 ~ cog_sim_pred_tics3_mn,
          time > 5 ~ cog_sim_pred_tics16_mn
        ),
        cog_blend_se = case_when(
          time <= 5 ~ cog_sim_pred_tics3_se,
          time > 5 ~ cog_sim_pred_tics16_se
        )
      ) %>% relocate(cog_blend_mn, .after = cog_sim_pred_tics16_mn)

r3taa <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics3_se:cog_blend_se,cog_sim_pred_tics16_mn)) %>% 
    pivot_longer(cols = c(cog_sim_pred_tics3_mn,cog_blend_mn),
        names_to = "label", values_to = "cog") %>% mutate(
          label = sub("_mn$","",label)
        )
r3tab <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics3_mn:cog_blend_mn,cog_sim_pred_tics16_se)) %>% 
    pivot_longer(cols = c(cog_sim_pred_tics3_se,cog_blend_se),
        names_to = "label", values_to = "cog_se") %>% mutate(
          label = sub("_se$","",label)
        )

r3ta <- r3taa %>% left_join((r3tab %>% 
    dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>% 
    mutate(
        cognitive_measure = factor(case_when(
            label == "cog_sim_pred_tics3" ~ "TICS (wave 3)",
            label == "cog_blend" ~ "Blended TICS (wave 3),TICS (wave 16)"
        ), levels = c("Blended TICS (wave 3),TICS (wave 16)","TICS (wave 3)"))
    )

r4ta <- r3ta %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
    mean_cog = mean(cog),
    se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
    lower = mean_cog - (1.96 * se_cog),
    upper = mean_cog + (1.96 * se_cog),
)

rm(r2,r3,r3a,r3b,r2ta,r3ta,r3taa,r3tab)


ggplot(data = r4, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    facet_wrap(vars(cognitive_measure)) +
    labs(caption = "True versus simulated cognitive trajectories (1000 Simulations) - All TICS Measures")

ggplot(data = r4ta, aes(x = time, y = mean_cog, color = slope_true_qrtl,
      linetype = cognitive_measure)) +
    geom_line() +
    geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
    labs(caption = str_wrap("TICS (wave 3) all assessments versus TICS (wave 3 - assessments bl-5) blended with TICS (wave 16 - assessments 6-10) (1000 Simulations)",width = 70))